自我代码提升之梯度提升树
作者：数据取经团——JQstyle
数据取经团（公众号：zlx19930503）
专注R、Python数据分析挖掘、可视化、机器学习
本文代码及数据获取方式：关注Python爱好者社区，在公众号中回复 GBDT
本期将为大家带来基于决策树基础模型的第二种算法，梯度提升树（即GBDT算法）。
GBDT的基本原理
梯度提升树属于Boosting集成学习算法的一种，其思想不同于随机森林、Bagging的并行化、投票的流程，GBDT模型所输出的结果是由其包含的若干棵决策树累加而得到的，每一棵子决策树都是实现对先前决策树组预测残差的拟合，是对先前模型的结果的一种“修正”。梯度提升树既可以用于回归问题（此时被称为CART回归树），也可以被用于解决分类问题（GBDT分类树），本文主要介绍用于分类的GBDT模型。
对于子模型类型的选择，GBDT方法和随机森林类似，通常均会选择CART回归树作为基模型。下面给出GBDT分类决策树一般的建模流程：
该流程即多分类GBDT算法，我们可以对其中各个步骤进行解读：
第0步：初始化各个类别下提升树的输出结果（可以将各个类别均设为0，也可以随机假设，需要注意的是，决策树输出的结果并不是分类结果，需要进行sigmoid函数转换），假定有K个类别，每个类别均有一个输出结果列，每条样本记录在各个列的对应取值为0或1，但对同一个样本，只能在一个类别上取1（该样本实际属于这个类别）其他均取0；
第1步：外层循环，所设定的决策树的数量为M，对其中任何一颗决策树进行循环；
第2步：对之前决策树在各个类别列上的输出结果进行函数转换，利用Softmax函数将这些结果转换为概率值，如初始化各列为0的情况下，各类别下的概率值进过Softmax转换以后均为1/K；
第3步：内层循环，对K个类别中的每一类进行循环；
第4步：求当前类别上残差的概率负梯度，Softmax所对应的损失函数los的负梯度即为y-p；
第5步：训练一棵包含J个叶节点的决策树（一般类似cart回归树），树的节点切分依据的指标为左子树样本和目标值的残差平方和+右子树样本和目标值的残差平方和-当前所有样本和目标值的残差平方和，最后叶节点输出当前样本的均值；
第6步：计算当前每个叶子节点的增益值，图中已经给出了公式
第7步：更新当前决策树组的输出公式，即之前的输出结果+当前决策树的输出结果
该流程将持续进行内外循环，直至结束
为了简化对问题的理解，本文所实现的GBDT是基于二分类问题的，因此本文将对上述流程进行细微的修改，即将多分类函数（通常为Softmax函数）改为适用于二分类的概率输出的Sigmoid函数。
对于二分类的GBDT的定义，通常会考虑将两个类别分别标记为1和-1进行处理，本文与之不同之处在于，将问题直接简化为0和1的分类问题，因此在对于损失函数的定义和梯度的求解会有所不同。另外，在GBDT的结果输出公式的定义上，对每个回归树的结果均会乘以一个步长step，使得延负梯度方向的前景不会过快而产生过拟合。
利用Sigmoid函数将GBDT决策树组的结果f(x)转化为一个概率值P(y=1)，该值越接近于1则说明越有可能被划分为1类，越接近于0则被归为0类（本人认为这样处理相比于-1和1的分类问题更便于理解，且梯度计算更为简单）。下面给出了本文所定义的损失函数和对应的梯度：
因此，最后本文给出的建模流程如下：
第0步：初始化提升树f(x)的输出结果直接设为0，此时sigmoid转化以后的概率是0.5；
第1步：循环，所设定的决策树的数量为M，对其中任何一颗决策树进行循环；
第2步：对之前决策树的输出结果进行sigmoid函数转换，将结果f(x)转换为概率值；
第3步：求当前类别上残差的概率负梯度；
第4步：对负梯度训练一棵回归树Rm，在训练过程中可以对树的深度depth、节点数J等进行控制；
第5步：计算当前每个叶子节点的增益值r；
第6步：更新当前决策树组的输出公式，即之前的输出结果+当前决策树的输出结果*step
第7步：重复2-6直至循环结束
GBDT的特点
相比于之前所介绍的诸多模型（KNN，逻辑回归等），GBDT属于较为复杂的预测模型之一，在很多现实的业务场景中已经得到运用。其优点在于：首先，基于cart树的GBDT算法可以很好的处理离散型和连续型数据指标，包括缺失值，并且可以运用于分类和回归问题；其次，作为一种集成学习的模型，GBDT往往具有更好的建模效果。相对的，GBDT的穿行算法流程在实践中效率往往较低，每一棵决策树的建立均需要依赖之前的结果，可并行化优化程度较低。
二分类GBDT的简单实现
接下来，我们尝试采用python实现GBDT算法，本次依旧采用马疝病horse数据集来进行测试，首先加相应的模块和数据集。代码见附件。
然后，我们定义一个cart回归树类cart_regression类，与之前的cart决策树不同的是，该模型每次节点分裂的依据是父节点样本均方误差与两个子节点均方误差和的差值，而预测的结果也是定量的节点样本均值。代码见附件。
接下来就开始定义GBDT类,该类是cart_regression的一个子类。在类中，首先定义sigmoid函数，用于将GBDT各个子模型的结果转化为0-1之间的概率。
class GBDT(cart_regression):
cart_list = []
def sigmoid(self,x):
return 1./(1+np.exp(-x))
然后重新定义cart子模型的训练函数，将其输出的结果定义成基于残差函数所及森的负梯度的增益的形式。
def cart_tree_fit(self,x,y,x_names,thre_num=5,
max_depth=3,now_depth=0): #重写CART树的训练函数
tree = {}
y_res = 1/2. * float(sum(y))/(sum(np.abs(y)*(1-np.abs(y))))
if y.shape[0] <= thre_num: #停止条件：当Y的数量低于阈值
return y_res
if now_depth >= max_depth:
return y_res
sel_x = 0
new_cut_x, break_point = self.new_cut_x(x,y) #提取各个X字段切分点和分享以后的X
gini_min = self.res_calc(x[:,0],y)
for j in range(1,new_cut_x.shape[1]): #找出误差最小的字段
if self.res_calc(new_cut_x[:,j],y) < gini_min:
gini_min = self.res_calc(new_cut_x[:,j],y)
sel_x = j
new_x_low = x[x[:,sel_x]<=break_point[sel_x],:]
new_x_high = x[x[:,sel_x]>break_point[sel_x],:]
new_y_low = y[x[:,sel_x]<=break_point[sel_x]] #当前字段低于切分点的Y数据集
new_y_high = y[x[:,sel_x]>break_point[sel_x]] #当前字段高于切分点的Y数据集
names_new = x_names#[np.array(range(x_names.shape[0]))!=sel_x]
label_low = x_names[sel_x] +'<=%s' %break_point[sel_x] #节点标签1
label_high = x_names[sel_x] +'>%s' %break_point[sel_x] #节点标签2
if np.unique(new_y_low).shape[0]<2 or np.unique(new_y_high).shape[0]<2:
return y_res
tree[label_low] = self.cart_tree_fit(new_x_low,new_y_low,names_new,thre_num,max_depth=max_depth,now_depth=now_depth+1) #子节点递归1
tree[label_high] = self.cart_tree_fit(new_x_high,new_y_high,names_new,thre_num,max_depth=max_depth,now_depth=now_depth+1) #子节点递归2
if tree[label_low] == tree[label_high]:
return tree[label_high]
return tree
与此同时，因为GBDT包含不止一个CART模型，因此我们需要对cart预测函数也进行重写，将模型对象也作为输入参数：
def cart_predict_line(self,x,x_names,model): #重写CART树的单行预测函数
import re
if isinstance(model,dict):
sel_x = x[x_names==re.split("<=|>",model.keys()[0])[0]]
bp = re.split("<=|>",model.keys()[0])[1]
if sel_x <= float(bp):
key_x = x_names[x_names==re.split("<=|>",model.keys()[0])[0]][0] + "<=" + bp
else:
key_x = x_names[x_names==re.split("<=|>",model.keys()[0])[0]][0] + ">" + bp
return self.cart_predict_line(x,x_names,model[key_x])
else:
return model
def cart_predict(self,x,x_names,model): #重写CART树预测函数
if model == {}:
return "Please fit the model"
else:
result = []
for i in range(x.shape[0]):
result.append(self.cart_predict_line(x[i,:],x_names,model))
result = np.array(result)
return result
下一步，我们定义GBDT初始化函数，将决策树的初始输出值定义为0
def init_value(self,y): return np.zeros(y.shape[0])
在定义GBDT的训练函数之前，因为GBDT的训练过程中需要不断地通过对先前回归树组的输出结果进行梯度和增益运算来更新模型，因此我们先定义预测函数（概率），该函数的功能包含了两部分，即在训练中用于不断输出当前回归树组对训练集的拟合结果，和在预测中对预测集的预测。在预测过程中，可以调节参数选择用于预测的树的数量（低于GBDT模型中子树的总数）：
def gbdt_predict(self,x,x_names,n_tree=-1,step=0.1):
y_init = self.init_value(x)
if self.cart_list == []:
return self.sigmoid(y_init)
elif n_tree <= 0:
for i in range(len(self.cart_list)):
pre_tmp = self.cart_predict(x,x_names,self.cart_list[i])
y_init += pre_tmp * step
return self.sigmoid(y_init)
else:
for i in range(n_tree):
pre_tmp = self.cart_predict(x,x_names,self.cart_list[i])
y_init += pre_tmp * step
return self.sigmoid(y_init)
之后定义GBDT训练的主函数，该函数需要传入训练数据集，字段名和01标签，默认每棵树最小分裂样本数为10，最大树深为5，子树的数量为3，步长为0.1，根据需要可以进行调整。
def GBDT_fit(self,x,y,x_names,thre_num=10,
max_depth=5,n_tree=3,step=0.1):
self.cart_list = []
y_pre = self.gbdt_predict(x,x_names,step=step)
for i in range(n_tree):
y_gd = y - y_pre
model_tmp = self.cart_tree_fit(x,y_gd,x_names,thre_num=thre_num,
max_depth=max_depth)
self.cart_list.append(model_tmp)
y_pre = self.gbdt_predict(x,x_names,n_tree=-1)
最后，由于先前的预测函数只输出概率值（即sigmoid的结果）,因此再定义类别预测函数，到此整个GBDT的类定义完成。
def gbdt_predict_type(self,x,x_names,n_tree=-1,thre_value=0.5,step=0.1):
y_pre = self.gbdt_predict(x,x_names,n_tree,step=step)
for i in range(y_pre.shape[0]):
if y_pre[i] > thre_value:
y_pre[i] = 1
else:
y_pre[i] = 0
return y_pre
接下来，我们对实现的GBDT算法进行测试，建立对马疝病生存的预测的模型。设定训练100棵子树，每课子树的最大深度为3，步长为0.1。在预测中，我们分别采用不同数量的子树来对训练集和测试集进行预测，观察预测的准确度：
#horse数据集测试
gbdt_model=GBDT() #实例化
gbdt_model.GBDT_fit(horse_train_x,horse_train_y,horse_name,thre_num=10,max_depth=3,n_tree=100) #模型训练
for j in [1,3,5,8,10,15,20,25,30,50,75,100]: #对比采用不同的数量的子树模型的预测的效果
gbdt_model_pre = gbdt_model.gbdt_predict_type(horse_test_x,horse_name,n_tree=j) #预测
gbdt_model_pre_train = gbdt_model.gbdt_predict_type(horse_train_x,horse_name,n_tree=j) #预测
ACC = sum(gbdt_model_pre == horse_test_y)/float(len(horse_test_y))
ACC_train = sum(gbdt_model_pre_train == horse_train_y)/float(len(horse_train_y))
print "the depth: %s, the n_tree: %s, the Accuracy to train: %s, the Accuracy to test: %s" % (3,j,ACC_train,ACC)
经过观察，我们发现，当深度仅为3的弱回归树学习器达到75棵时，在测试集的准确度可以达到82%，当前的测试效果且此时与训练集的效果基本持平。需要注意的是，在没有正则化等处理的情况下，一昧地增加回归树的数量或者深度容易陷入过拟合的问题。
基于Sklearn的GBDT建模
在末尾简单介绍一下运用机器学习模块sklearn简单的GBDT建模过程。sklearn的中的GradientBoostingClassifier集成学习GBDT方法包含了几个主要参数为：学习率（即步长），最小叶节点样本数，最大树深，子采样比例（每棵自回归树基于的训练样本抽样比例，用于避免过拟合），迭代次数等，读者可以自行尝试不同的参数组合：
#利用sklearn进行GBDT建模
from sklearn.ensemble import GradientBoostingClassifier #加载模块#实例化模型，传入参数（学习率：0.1，最小叶节点样本数：3，最小划分叶节点样本数：20，最大树深：7，子采样比例：0.8，迭代次数：40
GBDT_model = GradientBoostingClassifier(learning_rate=0.1, min_samples_split=3,
min_samples_leaf=20,max_depth=7, subsample=1,n_estimators=40)
GBDT_model.fit(horse_train_x,horse_train_y) #训练模型
y_pre_GBDT=GBDT_model.predict(horse_test_x) #预测模型
sum(y_pre_GBDT==horse_test_y)/float(horse_test_y.shape[0]) #准确度
拓展：浅谈XGBoost
XGBoost是目前在各类数据竞赛中被使用较为广泛的算法之一，已经被证明可以在诸多问题中取得非常好的预测效果。根据其开发者陈天琪在XGBoost论文中的叙述，XGBoost的基本思想差不多可以理解成以损失函数的二阶泰勒展开式作为其替代函数，求解其最小化（导数为0）来确定回归树的最佳切分点和叶节点输出数值（这一点和cart回归树不同）。此外，XGBoost通过在损失函数中引入子树数量和子树叶节点数值等，充分考虑到了正则化问题，能够有效避免过拟合。在效率上，XGBoost通过利用独特的近似回归树分叉点估计和子节点并行化等方式，加上二阶收敛的特性，建模效率较一般的GBDT有了大幅的提升。本文也附加了XGBoost建模的简单实践，有兴趣的读者可以去进行尝试：
#XGBoost
import xgboost as xgb
xgb_model = xgb.XGBClassifier(max_depth=2,n_estimators=100)
xgb_model.fit(horse_train_x,horse_train_y)
y_pre_xgb = xgb_model.predict(horse_test_x)
sum(y_pre_xgb==horse_test_y)/float(horse_test_y.shape[0])
最后，需要注意的是，本文所理解的二分类GBDT提升树算法是基于多分类简化而来的，而事实上二分类梯度提升树有着自己的损失函数和增益的计算方式（和本文存在一些差异），因此需要细致去理解的朋友可以参考网上相关的二分类梯度提升树的博客。
本文代码及数据获取方式：关注Python爱好者社区，在公众号中回复 GBDT
Python爱好者社区历史文章大合集：
Python爱好者社区历史文章列表（每周append更新一次）
关注后在公众号内回复“课程”即可获取：
小编的Python入门视频课程！！！
崔老师爬虫实战案例免费学习视频。
丘老师数据科学入门指导免费学习视频。
陈老师数据分析报告制作免费学习视频。
玩转大数据分析！Spark2.X+Python 精华实战课程免费学习视频。
丘老师Python网络爬虫实战免费学习视频。